# needed packages
pacman::p_load("brms", "glmmTMB", "metafor", "tidyverse",
"tidybayes", "orchaRd", "here")TUTORIAL: One dataset, four meta-analyses: synthesising mean effects, within-population variability, and between-population heterogeneity in ecology
Intended audience and scope of this tutorial
This tutorial is intended for researchers who already have a basic understanding of meta-analysis and some familiarity with statistical modelling in R. It is aimed at readers who wish to develop practical skills in location–scale meta-analytic modelling, using standard meta-analytic models as a foundation for understanding how mean effects and variability can be jointly analysed.
The scope of this tutorial is deliberately focused on conceptual comparison and practical illustration rather than exhaustive coverage. We do not provide an introduction to meta-analytic theory, effect size calculation, or general R programming, as these topics are well covered elsewhere. Instead, the tutorial uses a single dataset to demonstrate four complementary meta-analytic perspectives, highlighting their assumptions, strengths, and limitations.
Introduction
This online tutorial walks you through each step of meta-analysis and location-scale meta-regression using three R packages: metafor, brms, and glmmTMB. With one real dataset, we show how to analyse mean effects, within-population variability, and context-dependent heterogeneity, from data preparation through model interpretation.
Detailed methodological explanations and additional examples are beyond the scope of this tutorial - readers interested in further details may find the following website and relevant papers helpful:
Meta-analysis
Online tutorial: https://itchyshin.github.io/Meta-analysis_tutorial/
Related paper:
Nakagawa S, Yang Y, Macartney EL, Spake R, Lagisz M. (2023). Quantitative evidence synthesis: a practical guide on meta-analysis, meta-regression, and publication bias tests for environmental sciences. Environmental Evidence. 12:8. https://doi.org/10.1186/s13750-023-00301-6Location–scale meta-analysis models
Online tutorial: https://itchyshin.github.io/location-scale_meta-analysis/
Related paper:
Nakagawa S, Mizuno A, Morrison K, Ricolfi L, Williams C, Drobniak SM, Lagisz M, Yang Y. (2025). Location‐Scale Meta‐Analysis and Meta‐Regression as a Tool to Capture Large‐Scale Changes in Biological and Methodological Heterogeneity: A Spotlight on Heteroscedasticity. Global Change Biology. 31:e70204. https://doi.org/10.1111/gcb.70204Location–scale models
Online tutorial: https://ayumi-495.github.io/Eco_location-scale_model/
Related paper:
Nakagawa S, Ortega S, Gazzea E, Lagisz M, Lenz A, Lundgren E, Mizuno A. (2025). Location–scale models in ecology and evolution: Heteroscedasticity in continuous, count and proportion data. Methods in Ecology and Evolution. https://doi.org/10.1111/2041-210x.70203
Note
Outputs for the scale component in scale parts are reported on the log scale. Guidance on how to interpret these parameters, including back-transformation and interpretation is provided in the third online tutorial.
Setup
Load packages
First, we need to load necessary packages.
Important
Until recently. equalto(), which is required to conduct meta-analyses in glmmTMB, was not implemented. If you want to run meta-analytic models (including meta-regression), please ensure that you have installed the latest version of the package.
Also, we have created a custom function for visualising the outputs of meta-analysis models fitted in brms If you are interested, please expand the folded code below to view it.
Code
library(patchwork)
get_variables_dynamic <- function(model, patterns) {
variables <- get_variables(model)
idx <- unique(unlist(lapply(patterns, function(p) grep(p, variables))))
if (length(idx) == 0) return(character(0))
variables[idx]
}
# preparation
rename_vars <- function(variable) {
# fixed effects
variable <- gsub("b_Intercept", "b_l_intercept", variable)
variable <- gsub("b_sigma_Intercept", "b_s_intercept", variable)
variable <- gsub("b_Org.fertilizer.typeplant", "b_l_contrast", variable)
variable <- gsub("b_sigma_Org.fertilizer.typeplant", "b_s_contrast", variable)
# random effects
variable <- gsub("sd_study_ID__Intercept", "sd_study_ID", variable)
variable <- gsub("sigma", "sd_effect_ID", variable)
return(variable)
}
# Function to visualise fixed effects
visualize_fixed_effects <- function(model) {
fixed_effect_vars <- get_variables_dynamic(model, "^b_")
if (length(fixed_effect_vars) == 0) {
message("No fixed effects found")
return(NULL)
}
tryCatch({
fixed_effects_samples <- model %>%
spread_draws(!!!syms(fixed_effect_vars)) %>%
pivot_longer(cols = all_of(fixed_effect_vars), names_to = ".variable", values_to = ".value") %>%
mutate(.variable = rename_vars(.variable))
ggplot(fixed_effects_samples, aes(x = .value, y = .variable)) +
stat_halfeye(
normalize = "xy",
point_interval = "mean_qi",
fill = "lightcyan3",
color = "lightcyan4"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
labs(y = "Fixed effects", x = "Posterior values") +
theme_classic()
}, error = function(e) {
message("Error in visualize_fixed_effects: ", e$message)
return(NULL)
})
}
# Function to visualize random effects
visualize_random_effects <- function(model) {
random_effect_vars <- get_variables_dynamic(model, c("^sd_study", "^sigma"))
#random_effect_vars <- random_effect_vars[random_effect_vars != "sd_es_ID__Intercept"]
if (length(random_effect_vars) == 0) {
message("No random effects found")
return(NULL)
}
tryCatch({
random_effects_samples <- model %>%
spread_draws(!!!syms(random_effect_vars)) %>%
pivot_longer(cols = all_of(random_effect_vars), names_to = ".variable", values_to = ".value") %>%
mutate(.variable = rename_vars(.variable)) #%>%
#mutate(.value = .value) # leave SD as it is
ggplot(random_effects_samples, aes(x = .value, y = .variable)) +
stat_halfeye(
normalize = "xy",
point_interval = "mean_qi",
fill = "olivedrab3",
color = "olivedrab4"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
labs(y = "Random effects (SD)", x = "Posterior values") +
theme_classic()
}, error = function(e) {
message("Error in visualize_random_effects: ", e$message)
return(NULL)
})
}
# Function to visualize random effects
visualize_random_effects2 <- function(model) {
random_effect_vars <- get_variables_dynamic(model, c("^sd_study"))
#random_effect_vars <- random_effect_vars[random_effect_vars != "sd_es_ID__Intercept"]
if (length(random_effect_vars) == 0) {
message("No random effects found")
return(NULL)
}
tryCatch({
random_effects_samples <- model %>%
spread_draws(!!!syms(random_effect_vars)) %>%
pivot_longer(cols = all_of(random_effect_vars), names_to = ".variable", values_to = ".value") %>%
mutate(.variable = rename_vars(.variable)) #%>%
#mutate(.value = .value) # leave SD as it is
ggplot(random_effects_samples, aes(x = .value, y = .variable)) +
stat_halfeye(
normalize = "xy",
point_interval = "mean_qi",
fill = "olivedrab3",
color = "olivedrab4"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
labs(y = "Random effects (SD)", x = "Posterior values") +
theme_classic()
}, error = function(e) {
message("Error in visualize_random_effects: ", e$message)
return(NULL)
})
}All the outputs from the code chunks are saved in the Google drive folder (URL) to speed up the rendering process of this tutorial. If you want to run all the code chunks by yourself, please set eval: true for those code chunks with eval: false.
Data preparation
We use a dataset from Ponisio et al. (2015), which synthesises experimental studies comparing organic and conventional crop yields across a wide range of crop types, management practices, and geographic regions. The full meta-dataset comprises over 1,000 organic–conventional yield comparisons drawn from more than 100 primary studies; detailed metadata, inclusion criteria, and data descriptions are provided in the original publication and its associated supplementary materials.
A subset of the data is used here to facilitate computation and clarity of presentation while preserving the essential structure of the dataset
datfull <- read.csv(here::here("data", "ponisio2014dataset.csv"), header = TRUE)
# dim(datfull)
# [1] 1071 70 -> 1071 effect sizes and 70 columns
# head(datfull) # check column names
# create ID using Author + Year + Journal - gives unique identifier for each study
datfull$study_ID <- as.numeric(as.factor(paste(datfull$Author, datfull$Year, datfull$Journal, sep = "_")))
datfull$effect_size_ID <- 1:nrow(datfull)
# filter to the analysis subset
# Org.fertilizer.type = either plant or animal and also Crop.type = cereals
dat <- datfull %>%
dplyr::filter(Org.fertilizer.type %in% c("plant", "animal")) %>%
dplyr::filter(Crop.type == "cereals") %>%
droplevels()
# dim(dat) # 318 effect sizesEffect size calculation
Original study calculated and used log response ratio (lnRR) and its variance. We will also calculate log variance ratio (lnVR) and log coefficient of variation ratio (lnCVR), and their variances, respectively. metafor has a function escalc to calculate various effect sizes including lnVR and lnCVR.
dat <- escalc(measure = "VR",
m1i = Mean.Org,
sd1i = SD.Org,
n1i = N.Org,
m2i = Mean.conv,
sd2i = SD.conv,
n2i = N.Conv,
data = dat,
var.names = c("lnVR", "var.lnVR")
)
dat <- escalc(measure = "CVR",
m1i = Mean.Org,
sd1i = SD.Org,
n1i = N.Org,
m2i = Mean.conv,
sd2i = SD.conv,
n2i = N.Conv,
data = dat,
var.names = c("lnCVR", "var.lnCVR")
)Model 1: Meta-analysis of mean
We will first fit multilevel meta-analysis and meta-regression models to estimate the overall mean effect size using lnRR effect size. We will model effect sizes as the response variable, with random intercepts specified for study identity and effect-size identity to account for the hierarchical structure of the data and non-independence among observations.
metafor
When we want to run the multilevel meta-analysis using metafor, we can use the rma.mv function. Here, we specify lnRR as the response variable and var.lnRR as the sampling variance. Random effects are specified using the random argument, where we include random intercepts for both study_ID and effect_size_ID.
mod_lnRR <- rma.mv(lnRR,
var.lnRR,
random = list(~ 1 | study_ID,
~ 1 | effect_size_ID),
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnRR)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -109.8272 219.6545 225.6545 236.9312 225.7311
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0736 0.2712 36 no study_ID
# sigma^2.2 0.0830 0.2881 318 no effect_size_ID
#
# Test for Heterogeneity:
# Q(df = 317) = 15841.3280, p-val < .0001
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# -0.3204 0.0534 -5.9951 317 <.0001 -0.4255 -0.2152 ***
#
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
i2_ml(mod_lnRR)
# I2_Total I2_study_ID I2_effect_size_ID
# 98.11404 46.09489 52.01915 
The estimated average effect size was negative: \(\beta_0\) = -0.32, 95% confidence interval, CI = [-0.43, -0.22]), indicating that organic farming reduced crop yield by approximately 27% compared to conventional farming (1- exp(-0.32) = 1- 0.73). There was also high heterogeneity among effect sizes (\(I^2_{\text{Total}}\) = 98.11%), with contributions from both the study level and the effect size level (\(I^2_{\text{between}}\) = 46.09%, \(I^2_{\text{within}}\) = 52.02%). This indicates that effect sizes vary widely across studies as well as within studies.
Next, we fit a simple meta-regression to examine the type of fertiliser (animal- vs plant-based) which may explain variation in effect sizes.
mod_lnRR_mr <- rma.mv(lnRR,
var.lnRR,
mods = ~ 1 + Org.fertilizer.type,
random = list(~ 1 | study_ID,
~ 1 | effect_size_ID),
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnRR_mr)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -106.8735 213.7470 221.7470 236.7699 221.8756
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0614 0.2477 36 no study_ID
# sigma^2.2 0.0832 0.2885 318 no effect_size_ID
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 13892.9906, p-val < .0001
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 5.7785, p-val = 0.0168
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt -0.2663 0.0552 -4.8269 316 <.0001 -0.3749 -0.1578 ***
# Org.fertilizer.typeplant -0.1407 0.0585 -2.4038 316 0.0168 -0.2558 -0.0255 *
r2_ml(mod_lnRR_mr)
# R2_marginal R2_conditional
# 0.0323952 0.4429578 
The intercept represents the estimated mean lnRR for studies using animal-based fertilisers. The coefficient for plant-based fertilisers quantifies the difference in lnRR between plant- and animal-based fertilisers.
In the model summary, the estimated coefficient for plant-based fertilisers was negative and statistically significant, indicating that yields under plant-based organic fertilisers were, on average, lower than those under animal-based organic fertilisers. The moderator explained approximately 3.2% of the total variance in effect sizes (\(R^2_{\text{marginal}} = 0.032\)), indicating that differences between fertiliser types account for only a small proportion of the overall heterogeneity.
This pattern is illustrated in the orchard plot, which shows individual effect size estimates (lnRR) for plant- and animal-based fertilisers, with point sizes scaled by precision (1/SE). The black squares represent the model-predicted mean effects for each fertiliser type, with horizontal lines indicating their associated uncertainty intervals. Although the mean lnRR differs between fertiliser types, the substantial overlap and wide spread of individual effect sizes highlight that most of the variability remains unexplained by fertiliser type alone.
\(I^2\) vs \(R^2\)
\(I^2\) and \(R^2\) are both measures used in meta-analysis, but they quantify different aspects of variability in effect sizes.
\(I^2\) quantifies the proportion of the total observed variance in effect size estimates that is attributable to modelled heterogeneity rather than sampling variance. Its interpretation depends on the specified model structure and on which sources of heterogeneity are included, particularly in multilevel or location–scale meta-analytic models.
By contrast, \(R^2\) quantifies the proportion of variance in effect sizes that is explained by moderators in a meta-regression framework, conditional on the fitted model. \(R^2_{\text{marginal}}\) reflects the extent to which the included moderators explain heterogeneity, whereas \(R^2_{\text{conditional}}\) represents the proportion of variance explained by the full model, including both fixed effects and random effects, under the specified variance structure.
brms
As we mentioned in the main text, we can use the brms package to fit the same multilevel meta-analysis model. There are different ways to run these models, but in this example, we use the gr() function to set the sampling variance-covariance matrix. This approach is less common, but it can be much faster than usual (see more detailed information: link).
If you are not familiar with brms, we highly recommend checking out the instruction of brms first.
In the model below, sampling error is explicitly represented using a group-level term with a known variance–covariance structure. The term (1 | gr(effect_size_ID, cov = vcv)) specifies a random effect whose covariance matrix is fixed to the diagonal matrix of sampling variances for lnRR. This formulation is mathematically equivalent to supplying a sampling variance matrix in metafor::rma.mv(), but implemented here using brms.
Because sampling variances are assumed to be known in meta-analysis, the standard deviation of this group-level effect is fixed to 1. As a result, the magnitude of sampling error is fully determined by the supplied variance–covariance matrix and is not estimated from the data.
vcv <- diag(dat$var.lnRR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID
form_ma1 <- bf(lnRR ~ 1 +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
)
# generate default priors
prior_ma1 <- default_prior(form_ma1,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_ma1$prior[3] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1
prior_ma1
# fitting model
fit_ma1 <- brm(
formula = form_ma1,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_ma1,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_ma1)
# Family: gaussian
# Links: mu = identity
# Formula: lnRR ~ 1 + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.28 0.04 0.21 0.38 1.00 1839 3325
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept -0.32 0.06 -0.43 -0.21 1.00 1842 2814
#
# Further Distributional Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.29 0.01 0.26 0.32 1.00 7189 4944
#
# Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).
The estimated overall mean and between-study standard deviation is comparable to the study-level variance component obtained from the metafor model. In contrast, the residual parameter sigma represents unexplained heterogeneity at the effect-size level beyond sampling error, corresponding to the effect-size-level variance component in rma.mv.
One important point concerns how random effects are reported. In brms, random effect estimates are presented as standard deviations (SD), whereas metafor reports variances . Therefore, to compare estimates between the two packages, you need to square the standard deviation from brms to obtain the variance.
The estimated overall mean and between-study standard deviation are comparable to the study-level variance component obtained from the metafor model. In contrast, the residual parameter sigma represents unexplained heterogeneity at the effect-size level beyond sampling error, corresponding to the effect-size-level variance component in rma.mv().
One important point concerns how random effects are reported. In brms, random-effect estimates are presented as standard deviations (SD), whereas metafor reports variances. Therefore, to compare estimates between the two packages, the standard deviations reported by brms must be squared to obtain the corresponding variance estimates.
Finally, small numerical differences between estimates obtained from brms and metafor are expected, even when the underlying model structures are closely aligned. This is because brms fits models using Bayesian inference, which relies on probabilistic sampling to approximate posterior distributions, whereas metafor uses frequentist estimation. As a result, point estimates may differ slightly due to Monte Carlo error and differences in estimation frameworks, rather than reflecting substantive discrepancies between the models.
Let’s now fit a meta-regression model including fertiliser type as a moderator - the outputs are very similar to those from metafor.
vcv <- diag(dat$var.lnRR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID
form_mr1 <- bf(lnRR ~ 1 + Org.fertilizer.type +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
)
# generate default priors
prior_mr1 <- default_prior(form_mr1,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_mr1$prior[5] = "constant(1)" # meta-analysis assumes
# fitting model
fit_mr1 <- brm(
formula = form_mr1,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_mr1,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_mr1)
# Family: gaussian
# Links: mu = identity
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.26 0.04 0.19 0.35 1.00 2384 4173
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept -0.27 0.06 -0.38 -0.15 1.00 2787 3503
# Org.fertilizer.typeplant -0.14 0.06 -0.26 -0.02 1.00 4423 4398
#
# Further Distributional Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.29 0.01 0.26 0.32 1.00 6530 4799
#
# Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).You may wonder why we specify ~1 + Org.fertilizer.type rather than simply ~ Org.fertilizer.type. In fact, there is no difference between these two specifications. In R model formulas, an intercept term is included by default, so both formulas fit a model with an intercept.
We write 1 explicitly here for clarity. Doing so makes it clear that the model includes an overall intercept, which represents the baseline mean effect size for the reference level of the moderator. This is particularly helpful in a meta-analytic context, where the intercept has a clear interpretation as the overall mean effect. For this reason, we explicitly include 1 in all model formulas throughout this tutorial. This convention makes the role of the intercept transparent and helps avoid confusion when comparing different model specifications.
glmmTMB
glmmTMB is the popular package for fitting generalised linear mixed models. Unlike dedicated meta-analytic packages such as metafor, it does not automatically account for known sampling variances of effect sizes - so when we implement a meta-analysis in glmmTMB, we need explicitly specify the sampling error structure of the effect sizes.
In practice, this require three steps:
1. Constructing a variance-covariance matrix (VCV) that contains the known sampling variances of each effect size
vcv <- diag(dat$var.lnRR)
row.names(vcv) <- colnames(vcv) <- dat$effect_size_IDHere, we create a variance–covariance matrix whose diagonal elements are the known sampling variances of the effect sizes (var.lnRR). Off-diagonal elements are set to zero, reflecting the assumption that effect sizes are independent. This mirrors the standard assumption made in most meta-analyses.
2. Treating each effect size as a unique random-effect level
dat$effect_size_ID <- factor(dat$effect_size_ID)In glmmTMB, grouping variables for random effects must be factors. By converting effect_size_ID to a factor, we ensure that each effect size is treated as a distinct level in the random-effects structure. This step is required for the model to correctly associate each observation with its corresponding sampling variance.
3. Constraining the corresponding random-effect variance structure to be equal to the known VCV
g <- rep(1, nrow(dat))The equalto() term in glmmTMB requires a grouping factor, even when all observations belong to a single group. The dummy variable g serves this purpose. All observations are assigned to the same group, while the variance structure across effect sizes is fully defined by the supplied VCV matrix.
Conceptually, this step allows us to include sampling error as a random effect with a fixed, known variance structure.
You can also refer to this supplementary material for more details: https://coraliewilliams.github.io/equalto_sim_study/webpage.html.
glmm_lnRR_ma <- glmmTMB(lnRR ~ 1 +
(1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv), # incorporates the known sampling variances of each effect size
data = dat,
REML = TRUE
)
summary(glmm_lnRR_ma)
# Family: gaussian ( identity )
# Formula: lnRR ~ 1 + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 231.4 242.7 -112.7 225.4 315
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.0735694 0.27124
# g effect_size_ID33 0.0047213 0.06871
# effect_size_ID34 0.0426425 0.20650 0.00
# effect_size_ID35 0.0698313 0.26426 0.00 0.00
# effect_size_ID36 0.0798567 0.28259 0.00 0.00 0.00
# effect_size_ID37 0.0837435 0.28938 0.00 0.00 0.00 0.00
# effect_size_ID39 0.0095870 0.09791 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.0091520 0.09567 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.0363323 0.19061 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.0331452 0.18206 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.0012338 0.03513 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.0830247 0.28814
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.083
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.32039 0.05358 -5.98 2.23e-09 ***In short, you should interpret the output as follows:
Overall effect estimate:
(Intercept)Between-study variance estimate:
study_ID (Intercept)varianceEffect-size–level (sampling) variance:
Residual/sigma^2
All estimate is identical to metafor output, and also close to brms one.
NOTE:
The grouping factor g is a technical dummy variable required by glmmTMB when using equalto() to impose a known variance-covariance structure. All observations belong to a single level of g, and the variance structure across effect sizes is entirely defined by the supplied VCV matrix.
As a result, the multiple variance components listed under g in the model output correspond to the known sampling variances of individual effect sizes and are not estimated parameters. These components should not be interpreted. Only the study-level variance and the residual variance are of substantive interest.
When you want to run the meta-regression, you will simply add the moderators in the model:
glmm_lnRR_mr <- glmmTMB(lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
data = dat,
REML = TRUE
)
summary(glmm_lnRR_mr)
# Family: gaussian ( identity )
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +
# effect_size_ID | g, vcv)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 231.9 246.9 -111.9 223.9 314
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.0613541 0.24770
# g effect_size_ID33 0.0047213 0.06871
# effect_size_ID34 0.0426425 0.20650 0.00
# effect_size_ID35 0.0698313 0.26426 0.00 0.00
# effect_size_ID36 0.0798567 0.28259 0.00 0.00 0.00
# effect_size_ID37 0.0837435 0.28938 0.00 0.00 0.00 0.00
# effect_size_ID39 0.0095870 0.09791 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.0091520 0.09567 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.0363323 0.19061 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.0331452 0.18206 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.0012338 0.03513 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.0832439 0.28852
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.0832
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.26633 0.05522 -4.823 1.41e-06 ***
# Org.fertilizer.typeplant -0.14069 0.06008 -2.342 0.0192 * Great - the estimates match those from metafor and are also very close to the results from brms :)
Model 2: Meta-analysis of variability
In this section we analyse variability effects using lnCVR and lnVR.
lnCVR quantifies the log ratio of coefficients of variation (CVs), that is, changes in relative variability after accounting for mean differences.
lnVR quantifies the log ratio of standard deviations (SDs), that is, changes in absolute variability between groups.
metafor
lnCVR
mod_lnCVR <- rma.mv(lnCVR,
var.lnCVR,
random = list(~ 1 | study_ID,
~ 1 | effect_size_ID),
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnCVR)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -439.8288 879.6576 885.6576 896.9343 885.7342
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0602 0.2455 36 no study_ID
# sigma^2.2 0.2976 0.5455 318 no effect_size_ID
#
# Test for Heterogeneity:
# Q(df = 317) = 803.2915, p-val < .0001
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# 0.3430 0.0709 4.8354 317 <.0001 0.2034 0.4825 ***
#
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
i2_ml(mod_lnCVR)
# I2_Total I2_study_ID I2_effect_size_ID
# 56.525070 9.517436 47.007634
The intercept from the multilevel meta-analysis represents the overall mean lnCVR across all comparisons. A positive estimate indicates that the treatment group has a larger CV than the control group, meaning greater relative variability (lower stability) after accounting for potential differences in the mean values.
Here, the estimated lnCVR is \(\beta_0\) = 0.34 (95% CI= [0.20, 0.48]), suggesting that organic farming is associated with approximately 41% higher relative variability in crop yields compared to conventional farming (= exp(0.34)−1).
\(I^2\) values indicate that a substantial portion of the total heterogeneity (\(I^2_{\text{Total}}\) = 56.53%) arises from variability among effect sizes within studies (\(I^2_{\text{within}}\) = 47.01%), while a smaller share is due to variability between studies (\(I^2_{\text{between}}\) = 9.52%).
lnVR
mod_lnVR <- rma.mv(lnVR,
var.lnVR,
random = list(~ 1 | study_ID,
~ 1 | effect_size_ID),
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnVR)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -423.3326 846.6651 852.6651 863.9418 852.7418
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0640 0.2530 36 no study_ID
# sigma^2.2 0.1886 0.4342 318 no effect_size_ID
#
# Test for Heterogeneity:
# Q(df = 317) = 658.8664, p-val < .0001
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# 0.0198 0.0676 0.2924 317 0.7702 -0.1133 0.1528
#
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
i2_ml(mod_lnVR)
# I2_Total I2_study_ID I2_effect_size_ID
# 50.62070 12.82949 37.79121 
The overall lnVR estimate is close to zero and the 95% CI overlaps zero (\(\beta_0\) = 0.02, 95% CI = [-0.11, 0.15]), indicating little evidence for a systematic difference in SD between the two groups. On the original scale, this corresponds to an SD ratio close to one, suggesting that the magnitude of within-group variation is broadly similar between the two comparison groups..
Despite the absence of a clear average lnVR effect, heterogeneity is substantial. The test for heterogeneity is highly significant, and approximately half of the observed variability in lnVR cannot be attributed to sampling error alone (\(I^2 \approx\) 50.6%). Only a smaller proportion of this heterogeneity is attributable to between-study differences (\(I^2 \approx\) 12.8%), whereas a larger share arises at the effect-size level within studies (\(I^2 \approx\) 37.8%). This pattern indicates considerable context dependence or within-study variation in absolute variability effects, even though the average effect across all comparisons is close to zero.
When considered alongside lnCVR, a clearer picture emerges. While lnVR shows little evidence for a consistent difference in absolute variability, lnCVR indicates a positive average effect, suggesting higher relative variability (CV) in the treatment group. This divergence arises because lnCVR explicitly standardises variability by the mean, whereas lnVR does not. When group means differ, absolute variability may remain similar (lnVR \(\approx\) 0) while relative variability increases (lnCVR > 0).
In practical terms, this pattern implies that outcomes may scale with the mean. Although the absolute spread of values is comparable across groups, the spread is larger relative to the expected outcome under the treatment condition. Consequently, differences in variability are more apparent on the mean-standardised scale than on the absolute scale.
lnCVR
mod_lnCVR_mr <- rma.mv(lnCVR,
var.lnCVR,
mods = ~ Org.fertilizer.type,
random = list(~ 1 | study_ID,
~ 1 | effect_size_ID),
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnCVR_mr)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -438.7285 877.4570 885.4570 900.4800 885.5856
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0602 0.2454 36 no study_ID
# sigma^2.2 0.3001 0.5478 318 no effect_size_ID
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 803.2270, p-val < .0001
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 0.3026, p-val = 0.5826
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.3150 0.0874 3.6020 316 0.0004 0.1429 0.4870 ***
# Org.fertilizer.typeplant 0.0644 0.1171 0.5501 316 0.5826 -0.1660 0.2949
r2_ml(mod_lnCVR_mr)
# R2_marginal R2_conditional
# 0.00281066 0.16943708 
lnVR
mod_lnVR_mr <- rma.mv(lnVR,
var.lnVR,
mods = ~ Org.fertilizer.type,
random = list(~ 1 | study_ID,
~ 1 | effect_size_ID),
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnVR_mr)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -421.3055 842.6109 850.6109 865.6339 850.7396
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2.1 0.0560 0.2367 36 no study_ID
# sigma^2.2 0.1896 0.4355 318 no effect_size_ID
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 651.1955, p-val < .0001
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 2.4353, p-val = 0.1196
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.0876 0.0791 1.1075 316 0.2689 -0.0680 0.2433
# Org.fertilizer.typeplant -0.1644 0.1054 -1.5605 316 0.1196 -0.3717 0.0429
r2_ml(mod_lnVR_mr)
# R2_marginal R2_conditional
# 0.02621268 0.24827290 
Including fertiliser type as a moderator provided little evidence that it systematically explains variability effects on either scale. Here, the moderator coefficient represents the difference between plant-based and animal-based fertilisers (plant − animal).
For lnCVR, the estimated moderator effect was small and highly uncertain (\(\beta_{\text{plant-animal}}\) = 0.064, 95% CI = [-0.17, 0.29]), indicating that relative variability does not differ consistently between fertiliser categories.
Similarly, for lnVR, there was no clear support for fertiliser-type differences in absolute variability (\(\beta_{\text{plant-animal}}\) = −0.164, 95% CI = [-0.37, 0.04]). Although the point estimate suggested a tendency toward lower lnVR for plant-based fertilisers relative to animal-based fertilisers, this effect was highly uncertain and not statistically supported.
In both models, substantial residual heterogeneity remained after accounting for fertiliser type, with large effect-size–level variance components (lnCVR: \(\sigma^2\) = 0.30; lnVR: \(\sigma^2\) = 0.19). This indicates that neither relative nor absolute variability effects are well captured by this single moderator.
brms
When fitted with equivalent model structures, the metafor and brms implementations yield nearly identical estimates for both the overall effects and variance components.
lnCVR
vcv <- diag(dat$var.lnCVR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID
form_ma2 <- bf(lnCVR ~ 1 +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
)
# generate default priors
prior_ma2 <- default_prior(form_ma2,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_ma2$prior[3] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1
# fitting model
fit_ma2 <- brm(
formula = form_ma2,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_ma2,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_ma2)
# Family: gaussian
# Links: mu = identity
# Formula: lnCVR ~ 1 + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.27 0.08 0.12 0.45 1.00 1404 2356
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 0.34 0.07 0.20 0.49 1.00 3888 4078
#
# Further Distributional Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.55 0.05 0.45 0.66 1.00 896 1820
lnVR
vcv <- diag(dat$var.lnVR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID
form_ma3 <- bf(lnVR ~ 1 +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
)
# generate default priors
prior_ma3 <- default_prior(form_ma3,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_ma3$prior[3] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1
# fitting model
fit_ma3 <- brm(
formula = form_ma3,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_ma3,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_ma3)
# Family: gaussian
# Links: mu = identity
# Formula: lnVR ~ 1 + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.27 0.07 0.15 0.42 1.00 1458 2196
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 0.02 0.07 -0.12 0.17 1.00 2816 3765
#
# Further Distributional Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.44 0.06 0.32 0.56 1.01 469 856
When specified with equivalent fixed and random effects and the same sampling variance structure, the brms and metafor meta-regression models yield nearly identical results. Differences between outputs reflect estimation framework (Bayesian vs frequentist) rather than substantive differences in model inference.
lnCVR
vcv <- diag(dat$var.lnCVR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID
form_mr2 <- bf(lnCVR ~ 1 + Org.fertilizer.type +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
)
# generate default priors
prior_mr2 <- default_prior(form_mr2,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_mr2$prior[5] = "constant(1)" # meta-analysis assumes
# fitting model
fit_mr2 <- brm(
formula = form_mr2,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_mr2,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_mr2)
# Family: gaussian
# Links: mu = identity
# Formula: lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.27 0.08 0.13 0.43 1.00 1515 2605
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 0.32 0.09 0.14 0.50 1.00 4506 4337
# Org.fertilizer.typeplant 0.06 0.12 -0.18 0.30 1.00 3709 4471
#
# Further Distributional Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.55 0.05 0.45 0.66 1.00 987 1853lnVR
vcv <- diag(dat$var.lnVR)
rownames(vcv) <- colnames(vcv) <- dat$effect_size_ID
form_mr3 <- bf(lnVR ~ 1 + Org.fertilizer.type +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)) # this is m (sampling error)
)
# generate default priors
prior_mr3 <- default_prior(form_mr3,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_mr3$prior[5] = "constant(1)" # meta-analysis assumes
# fitting model
fit_mr3 <- brm(
formula = form_mr3,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_mr3,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_mr3)
# Family: gaussian
# Links: mu = identity
# Formula: lnVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.24 0.07 0.11 0.39 1.00 1038 1472
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 0.09 0.08 -0.06 0.25 1.00 3110 4251
# Org.fertilizer.typeplant -0.17 0.11 -0.38 0.04 1.00 2781 3712
#
# Further Distributional Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.44 0.06 0.32 0.56 1.00 443 659glmmTMB
We can also fit the same meta-analytic and meta-regression models using glmmTMB. The outputs are very similar to those from metafor and brms.
First, we need to prepare the variance-covariance matrix for the sampling errors as shown below…
# lnCVR
vcv <- diag(dat$var.lnCVR)
row.names(vcv) <- colnames(vcv) <- dat$effect_size_ID
dat$effect_size_ID <- factor(dat$effect_size_ID)
g <- rep(1, nrow(dat))
# lnVR
vcv <- diag(dat$var.lnVR)
row.names(vcv) <- colnames(vcv) <- dat$effect_size_ID
dat$effect_size_ID <- factor(dat$effect_size_ID)
g <- rep(1, nrow(dat))glmm_lnCVR <- glmmTMB(lnCVR ~ 1 + (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
data = dat,
REML = TRUE
)
summary(glmm_lnCVR)
# Family: gaussian ( identity )
# Formula: lnCVR ~ 1 + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 892.6 903.8 -443.3 886.6 315
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.06173 0.2485
# g effect_size_ID33 0.11111 0.3333
# effect_size_ID34 0.25000 0.5000 0.00
# effect_size_ID35 0.25000 0.5000 0.00 0.00
# effect_size_ID36 0.25000 0.5000 0.00 0.00 0.00
# effect_size_ID37 0.25000 0.5000 0.00 0.00 0.00 0.00
# effect_size_ID39 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.31923 0.5650
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.319
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.34418 0.07145 4.817 1.46e-06 ***glmm_lnVR <- glmmTMB(lnVR ~ 1 + (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
data = dat,
REML = TRUE
)
summary(glmm_lnVR)
# Family: gaussian ( identity )
# Formula: lnVR ~ 1 + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 858.4 869.7 -426.2 852.4 315
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.06401 0.2530
# g effect_size_ID33 0.11111 0.3333
# effect_size_ID34 0.25000 0.5000 0.00
# effect_size_ID35 0.25000 0.5000 0.00 0.00
# effect_size_ID36 0.25000 0.5000 0.00 0.00 0.00
# effect_size_ID37 0.25000 0.5000 0.00 0.00 0.00 0.00
# effect_size_ID39 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.18856 0.4342
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.189
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.01977 0.06787 0.291 0.771glmm_lnCVR_mr <- glmmTMB(lnCVR ~ 1
+ Org.fertilizer.type
+ (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
data = dat,
REML = TRUE
)
summary(glmm_lnCVR_mr)
# Family: gaussian ( identity )
# Formula: lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 896.6 911.7 -444.3 888.6 314
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.06125 0.2475
# g effect_size_ID33 0.11111 0.3333
# effect_size_ID34 0.25000 0.5000 0.00
# effect_size_ID35 0.25000 0.5000 0.00 0.00
# effect_size_ID36 0.25000 0.5000 0.00 0.00 0.00
# effect_size_ID37 0.25000 0.5000 0.00 0.00 0.00 0.00
# effect_size_ID39 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.32171 0.5672
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.322
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.31394 0.08819 3.560 0.000371 ***
# Org.fertilizer.typeplant 0.07016 0.11845 0.592 0.553664 glmm_lnVR_mr <- glmmTMB(lnVR ~ 1
+ Org.fertilizer.type
+ (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
data = dat,
REML = TRUE
)
summary(glmm_lnVR_mr)
# Family: gaussian ( identity )
# Formula: lnVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 860.7 875.8 -426.4 852.7 314
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.05602 0.2367
# g effect_size_ID33 0.11111 0.3333
# effect_size_ID34 0.25000 0.5000 0.00
# effect_size_ID35 0.25000 0.5000 0.00 0.00
# effect_size_ID36 0.25000 0.5000 0.00 0.00 0.00
# effect_size_ID37 0.25000 0.5000 0.00 0.00 0.00 0.00
# effect_size_ID39 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.11111 0.3333 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.25000 0.5000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual 0.18965 0.4355
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion estimate for gaussian family (sigma^2): 0.19
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.08762 0.07921 1.106 0.269
# Org.fertilizer.typeplant -0.16443 0.10620 -1.548 0.122Models 3 and 4: Location-scale meta-regression models
We now extend the above meta-analyses to location–scale meta-regression models, which allow us to examine whether both the mean effect size and the variability of effect sizes depend on a moderator variable. Here, we use Org.fertilizer.type (plant-based vs animal-based organic fertiliser) as a moderator.
The location component answers the question:
Does the mean effect size (i.e. lnRR, lnCVR, lnVR) differ between plant-based and animal-based organic fertilisers?The scale component answers the question:
Does the residual heterogeneity (i.e. variability of effect sizes after accounting for random effects) differ between plant-based and animal-based organic fertilisers?
metafor
At present, metafor does not support fully multilevel location–scale meta-analytic models. However, it is possible to fit a simpler location–scale model using the scale argument in rma(). This approach allows the residual heterogeneity to vary as a function of moderators, but does not accommodate multilevel random-effects structures.
Below, we demonstrate two types of models:
Simple location–scale meta-regression
Multilevel heteroscedastic meta-regression
Note
It is best to treat the metafor models presented above as approximations to full multilevel location-scale meta-analytic models. Although these models are extremely useful and sometimes sufficient in practice, ecological datasets often show more complex hierarchical structures, such as multiple effect sizes per study, evolutionary independence, or other sources of non-independence. To explicitly consider such complexity, more flexible modelling frameworks are required. In the following two sections, we introduce location-scale meta-analytic models fitted using brms and glmmTMB. These approach arrow the specification of rich multilevel structures, and in particular, brms enables the fitting of fully Bayesian location-scale models with complex random effects and covariance structures (though if we want to run this type of model, we need a lot of data…)
The formulation using rma() provides a straightforward way to test whether variability in effect sizes depends on the moderator, but it treats all effect sizes as independent and does not include a multilevel structure. scale argument allows us to specify a model for the residual variance, enabling us to test whether the variability of effect sizes differs across levels of the moderator.
lnRR
mod_lnRR_ls <- rma(yi = lnRR,
vi = var.lnRR,
mods = ~ Org.fertilizer.type,
scale = ~ Org.fertilizer.type,
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnRR_ls)
# Location-Scale Model (k = 318; tau^2 estimator: REML)
#
# logLik deviance AIC BIC AICc
# -161.2053 322.4106 330.4106 345.4336 330.5393
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 13892.9906, p-val < .0001
#
# Test of Location Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 67.0786, p-val < .0001
#
# Model Results (Location):
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt -0.1796 0.0290 -6.1935 316 <.0001 -0.2366 -0.1225 ***
# Org.fertilizer.typeplant -0.3666 0.0448 -8.1902 316 <.0001 -0.4547 -0.2785 ***
#
# Test of Scale Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 19.9223, p-val < .0001
#
# Model Results (Scale):
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt -2.4057 0.1423 -16.9060 316 <.0001 -2.6856 -2.1257 ***
# Org.fertilizer.typeplant 0.8071 0.1808 4.4634 316 <.0001 0.4513 1.1629 *** lnCVR
mod_lnCVR_ls <- rma(yi = lnCVR,
vi = var.lnCVR,
mods = ~ Org.fertilizer.type,
scale = ~ Org.fertilizer.type,
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnCVR_ls)
# Location-Scale Model (k = 318; tau^2 estimator: REML)
#
# logLik deviance AIC BIC AICc
# -439.5156 879.0311 887.0311 902.0541 887.1597
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 803.2270, p-val < .0001
#
# Test of Location Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 1.0750, p-val = 0.3006
#
# Model Results (Location):
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.2401 0.0620 3.8719 316 0.0001 0.1181 0.3621 ***
# Org.fertilizer.typeplant 0.1026 0.0990 1.0368 316 0.3006 -0.0921 0.2974
#
# Test of Scale Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 10.1534, p-val = 0.0016
#
# Model Results (Scale):
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt -1.4908 0.2473 -6.0286 316 <.0001 -1.9773 -1.0043 ***
# Org.fertilizer.typeplant 1.0399 0.3263 3.1864 316 0.0016 0.3978 1.6820 ** lnVR
mod_lnVR_ls <- rma(yi = lnVR,
vi = var.lnVR,
mods = ~ Org.fertilizer.type,
scale = ~ Org.fertilizer.type,
data = dat,
test = "t",
method = "REML"
)
summary(mod_lnVR_ls)
# Location-Scale Model (k = 318; tau^2 estimator: REML)
#
# logLik deviance AIC BIC AICc
# -415.2505 830.5010 838.5010 853.5240 838.6296
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 651.1955, p-val < .0001
#
# Test of Location Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 4.7895, p-val = 0.0294
#
# Model Results (Location):
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.0417 0.0457 0.9106 316 0.3632 -0.0483 0.1317
# Org.fertilizer.typeplant -0.1906 0.0871 -2.1885 316 0.0294 -0.3620 -0.0193 *
#
# Test of Scale Coefficients (coefficient 2):
# F(df1 = 1, df2 = 316) = 11.9043, p-val = 0.0006
#
# Model Results (Scale):
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt -2.8652 0.6270 -4.5699 316 <.0001 -4.0987 -1.6316 ***
# Org.fertilizer.typeplant 2.3044 0.6679 3.4503 316 0.0006 0.9903 3.6185 *** The models suggest that plant-based fertilisers are associated with both lower mean effects and substantially greater variability across studies. In particular, the scale component indicates that effect sizes are markedly more dispersed for plant-based fertilisers across all three variability metrics. However, these models treat all effect sizes as independent and do not account for the multilevel structure of the data.
Heteroscedastic meta-regression models allow residual variance to differ across levels of a moderator while incorporating a multilevel random-effects structure. By relaxing the assumption of constant heterogeneity, these models permit different levels of a moderator to be associated with different amounts of residual variability. Although this approach does not constitute a full multilevel location–scale model, it provides a useful intermediate step toward more flexible variance modelling. In particular, it enables the examination of moderator effects on both the mean effect size and the residual heterogeneity while accounting for the multilevel structure of the data.
In the present heteroscedastic meta-regression, the random-effects structure ~ Org.fertilizer.type | effect_size_ID allows residual heterogeneity at the effect-size level to differ across fertiliser types. By specifying struct = "DIAG", the model estimates fertiliser-type-specific variances without assuming covariance between categories, thereby relaxing the assumption of homogeneous residual variability while retaining a parsimonious specification.
Tip
We used this heteroscedastic meta-regression models for figures 3A and 3B in the main text.
lnRR
mod_lnRR_mr2 <- rma.mv(lnRR,
var.lnRR,
mods = ~ Org.fertilizer.type,
random = list(~ 1 | study_ID,
~ Org.fertilizer.type | effect_size_ID),
data = dat,
test = "t",
method = "REML",
struct = 'DIAG'
)
summary(mod_lnRR_mr2)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -105.4125 210.8250 220.8250 239.6037 221.0185
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2 0.0608 0.2466 36 no study_ID
#
# outer factor: effect_size_ID (nlvls = 318)
# inner factor: Org.fertilizer.type (nlvls = 2)
#
# estim sqrt k.lvl fixed level
# tau^2.1 0.0666 0.2581 134 no animal
# tau^2.2 0.0937 0.3061 184 no plant
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 13892.9906, p-val < .0001
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 5.9453, p-val = 0.0153
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt -0.2679 0.0537 -4.9852 316 <.0001 -0.3736 -0.1622 ***
# Org.fertilizer.typeplant -0.1399 0.0574 -2.4383 316 0.0153 -0.2528 -0.0270 * 
lnCVR
mod_lnCVR_mr2 <- rma.mv(lnCVR,
var.lnCVR,
mods = ~ Org.fertilizer.type,
random = list(~ 1 | study_ID,
~ Org.fertilizer.type | effect_size_ID),
data = dat,
test = "t",
method = "REML",
struct = 'DIAG'
)
summary(mod_lnCVR_mr2)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -434.6673 869.3347 879.3347 898.1134 879.5282
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2 0.0578 0.2405 36 no study_ID
#
# outer factor: effect_size_ID (nlvls = 318)
# inner factor: Org.fertilizer.type (nlvls = 2)
#
# estim sqrt k.lvl fixed level
# tau^2.1 0.1915 0.4376 134 no animal
# tau^2.2 0.5419 0.7361 184 no plant
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 803.2270, p-val < .0001
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 0.4957, p-val = 0.4819
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.3021 0.0821 3.6816 316 0.0003 0.1406 0.4635 ***
# Org.fertilizer.typeplant 0.0845 0.1200 0.7040 316 0.4819 -0.1517 0.3207 
lnVR
mod_lnVR_mr2 <- rma.mv(lnVR,
var.lnVR,
mods = ~ Org.fertilizer.type,
random = list(~ 1 | study_ID,
~ Org.fertilizer.type | effect_size_ID),
data = dat,
test = "t",
method = "REML",
struct = 'DIAG'
)
summary(mod_lnVR_mr2)
# Multivariate Meta-Analysis Model (k = 318; method: REML)
#
# logLik Deviance AIC BIC AICc
# -409.8610 819.7221 829.7221 848.5008 829.9156
#
# Variance Components:
#
# estim sqrt nlvls fixed factor
# sigma^2 0.0588 0.2425 36 no study_ID
#
# outer factor: effect_size_ID (nlvls = 318)
# inner factor: Org.fertilizer.type (nlvls = 2)
#
# estim sqrt k.lvl fixed level
# tau^2.1 0.0291 0.1707 134 no animal
# tau^2.2 0.5092 0.7136 184 no plant
#
# Test for Residual Heterogeneity:
# QE(df = 316) = 651.1955, p-val < .0001
#
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 316) = 1.9172, p-val = 0.1671
#
# Model Results:
#
# estimate se tval df pval ci.lb ci.ub
# intrcpt 0.0666 0.0713 0.9329 316 0.3516 -0.0738 0.2069
# Org.fertilizer.typeplant -0.1526 0.1102 -1.3846 316 0.1671 -0.3695 0.0643 
In the heteroscedastic meta-regression for lnRR, residual heterogeneity at the effect-size level differed between fertiliser types, with larger heterogeneity for plant-based fertilisers (\(\tau^2\) = 0.094) than for animal-based fertilisers (\(\tau^2\) = 0.067). The location component indicated a significantly lower mean lnRR for plant-based fertilisers relative to animal-based fertilisers (estimate = −0.14, 95% CI = [-0.25, -0.03]). Thus, even after accounting for multilevel structure, plant-based fertilisers were associated with both lower mean effects and moderately greater residual heterogeneity.
For lnCVR, allowing residual heterogeneity to vary across fertiliser types revealed substantial differences in heterogeneity magnitude. Effect-size–level heterogeneity was markedly larger for plant-based fertilisers (\(\tau^2\) = 0.54) than for animal-based fertilisers (\(\tau^2\) = 0.19). However, the location effect of fertiliser type on mean lnCVR remained small and statistically unsupported (estimate = 0.08, 95% CI = [-0.15, 0.32]), indicating no consistent difference in average relative variability between fertiliser categories.
A similar pattern was observed for lnVR. Residual heterogeneity was considerably larger for plant-based fertilisers (\(\tau^2\) = 0.51) than for animal-based fertilisers (\(\tau^2\) = 0.03), suggesting pronounced differences in the dispersion of absolute variability across fertiliser types. In contrast, the location effect on mean lnVR was not statistically supported (estimate = −0.15, 95% CI = [-0.37, 0.06]), indicating weak evidence for systematic differences in average absolute variability.
brms
brms provides a flexible framework for fitting location–scale (LS) models, in which the mean structure of the response variable (the location component) and the residual variability around that mean (the scale component) are modelled simultaneously. This allows us to examine not only how predictors affect average effect sizes, but also how they influence heterogeneity.
In this tutorial, we focus primarily on a location–scale meta-regression model with fixed effects in the scale component. This model represents the most stable and interpretable formulation and serves as the main working example throughout the brms section.
For completeness, we also briefly introduce two more complex extensions:
Location–scale meta-regression with fixed effects in the scale component.
This is the primary model used in this tutorial and the only one for which we provide full interpretation of results.Location–scale meta-regression with random effects in the scale component.
This model allows heterogeneity to vary across higher-level grouping factors (for example, studies), but may be difficult to fit reliably in practice.Location–scale meta-regression with correlated random effects in the location and scale components.
This formulation allows the mean and heterogeneity to be evolutionarily or structurally linked, but often suffers from convergence issues in applied meta-analytic datasets.
Models B and C are included to illustrate the range of model structures that are theoretically possible in brms. We provide one concrete example of each, but we do not interpret their results in detail, as these models frequently exhibit convergence problems and require a lot of data. If you are interested in more detail, please check Nakagawa et al. (2025a, b).
It is worth noting that, in the models presented here, we specify the scale component using the sigma ~ syntax in brms. This choice reflects the use of a Gaussian response distribution, for which sigma represents the residual standard deviation. For other response distributions, the scale component may be parametrised differently in brms (for example, using shape or phi - see here), and the corresponding model formula should be adapted accordingly.
Practical guidance
Models with random effects in the scale component, and especially models with correlated random effects across location and scale, require:
Sufficient numbers of studies
Adequate replication within studies
Careful prior specification
Thorough convergence diagnostics
For these reasons, we treat Models A and B as examples. We recommend that readers begin with a simpler location–scale model that uses fixed effects in the scale part, and only add complexity if the data clearly support it, and also, you have enough data.
In this model, we include fixed effects in both location and scale components, but no random effects in the scale component (of course, we can include any random effects on location component).
Between-study heterogeneity in mean effect sizes is captured through a random intercept for study_ID in the location component, and residual heterogeneity is allowed to vary systematically with predictors via fixed effects in the scale component. No random variation in residual variance across studies is assumed.
lnRR
form_ls_ma1 <- bf(lnRR ~ 1 + Org.fertilizer.type +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)), # this is m (sampling error)
sigma ~ 1 + Org.fertilizer.type # fixed effects only in scale component
)
# generate default prior
prior_ls_ma1 <- default_prior(form_ls_ma1,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_ls_ma1$prior[5] = "constant(1)" # fixing the variance to 1
# fitting model
fit_ls_ma1 <- brm(
formula = form_ls_ma1,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_ls_ma1,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_ls_ma1)
# Family: gaussian
# Links: mu = identity; sigma = log
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# sigma ~ 1 + Org.fertilizer.type
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.45 0.27 0.02 1.02 1.00 865 2021
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.26 0.04 0.19 0.35 1.00 2342 3505
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept -0.26 0.06 -0.38 -0.16 1.00 2437 3607
# sigma_Intercept -1.30 0.07 -1.45 -1.15 1.00 3773 3686
# Org.fertilizer.typeplant -0.15 0.06 -0.27 -0.03 1.00 5617 4995
# sigma_Org.fertilizer.typeplant 0.17 0.09 -0.01 0.35 1.00 9838 4318For lnRR, the location component indicated a negative average effect under animal-based fertilisers (\(\beta^{(l)}_{\text{animal}}\) = -0.26, 95% credible interval CrI = [−0.38, −0.16]), implying that organic yields are, on average, lower than conventional yields under animal-based fertilisation. The contrast coefficient for fertiliser type was also negative (\(\beta^{(l)}_{\text{animal-plant}}\) = -0.15, 95% CrI [−0.27, −0.03]), indicating an additional yield penalty under plant-based fertilisers relative to animal-based fertilisers. On the response-ratio scale, these estimates correspond to an expected organic yield of \(\exp(\beta_0)\) = \(\exp\)(-0.26) \(\approx\) 0.77, or about 23% lower than conventional yields, under animal-based fertilisers (the reference category). For plant-based fertilisers, the predicted response ratio is given by \(\exp(\beta_0\) + \(\beta_{\text{plant}}\)) = \(\exp\)(-0.26 - 0.15) \(\approx 0.66\), corresponding to organic yields that are approximately 34% lower than conventional yields.
The scale component, modelled on the log-\(\sigma\) scale, suggested a tendency towards greater residual heterogeneity under plant-based fertilisers (\(\beta^{(s)}_{\text{animal-plant}}\) = 0.17, 95% CrI = [−0.01, 0.35]). This corresponds to an approximate 19% increase in residual heterogeneity for plant-based systems ((0.17)-1 \(\approx\) 0.19), although the credible interval still overlapped zero, indicating moderate uncertainty in the magnitude of this effect.
lnCVR
form_ls_ma2 <- bf(lnCVR ~ 1 + Org.fertilizer.type +
(1|study_ID) + # this is u_l (the between-study effect)
(1|gr(effect_size_ID, cov = vcv)), # this is m (sampling error)
sigma ~ 1 + Org.fertilizer.type
)
# generate default priors
prior_ls_ma2 <- default_prior(form_ls_ma2,
data=dat,
data2=list(vcv=vcv),
family=gaussian())
prior_ls_ma2$prior[5] = "constant(1)" # meta-analysis assumes
# fitting model
fit_ls_ma2 <- brm(
formula = form_ls_ma2,
data = dat,
data2 = list(vcv=vcv),
chains = 2,
cores = 2,
iter = 6000,
warmup = 3000,
prior = prior_ls_ma2,
control = list(adapt_delta=0.95, max_treedepth=15)
)
summary(fit_ls_ma2)
# Family: gaussian
# Links: mu = identity; sigma = log
# Formula: lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + (1 | gr(effect_size_ID, cov = vcv))
# sigma ~ 1 + Org.fertilizer.type
# Data: dat (Number of observations: 318)
# Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
# total post-warmup draws = 6000
#
# Multilevel Hyperparameters:
# ~effect_size_ID (Number of levels: 318)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 1.00 0.00 1.00 1.00 NA NA NA
#
# ~study_ID (Number of levels: 36)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.26 0.08 0.12 0.43 1.00 1516 2572
#
# Regression Coefficients:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 0.30 0.09 0.13 0.48 1.00 4592 4409
# sigma_Intercept -0.83 0.13 -1.09 -0.58 1.00 1518 2937
# Org.fertilizer.typeplant 0.08 0.13 -0.17 0.32 1.00 4734 4478
# sigma_Org.fertilizer.typeplant 0.51 0.18 0.15 0.87 1.00 1459 2679
#
# Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).
For lnCVR, the location component indicated that organic systems tend to exhibit higher mean-corrected variability under animal-based fertilisers (\(\beta^{(l)}_{\text{animal}}\) = 0.30, 95% CrI = [0.13, 0.48]), corresponding to an average CV ratio of (0.30) \(\approx\) 1.35 - roughly 35% higher variability in organic relative to conventional systems.
The contrast coefficient for fertiliser type in the location component was small and highly uncertain (\(\beta^{(l)}_{\text{animal-plant}}\) = 0.08, 95% CrI [−0.17, 0.32]), suggesting that fertiliser category does not systematically shift the average stability difference between organic and conventional systems.
In contrast, the scale component revealed clear evidence that fertiliser type strongly governs the predictability of lnCVR effects. Residual heterogeneity was substantially larger under plant-based fertilisers \((\beta^{(s)}_{\text{animal-plant}}\) = 0.51, 95% CrI [0.15, 0.87]), corresponding to a multiplicative increase in residual variability of approximately 68% ((0.51) \(\approx\) 1.68). Under this model, the implied residual standard deviation for lnCVR under plant-based fertilisers is \(\sigma \approx \exp\)(-0.83 + 0.51) \(\approx\) 0.73.
Thus, while the expected magnitude of the stability gap (lnCVR) is broadly similar across fertiliser categories, the evidence is far less transferable under plant-based fertilisers. Some studies show near parity in stability between organic and conventional systems, whereas others show much larger differences, indicating strong context dependence.
The models presented here extend the basic location–scale meta-regression by allowing study-level heterogeneity to operate not only on the mean effect size, but also on the residual variance.
These models are included for completeness and to illustrate the flexibility of brms. In practice, they are often difficult to fit reliably and should be treated as advanced extensions rather than default choices. Accordingly, we provide example code but do not interpret the results in detail.
For transparency, examples in which the models were run but failed to converge are also documented in this online tutorial. You can also find more examples of complex location-scale models in another online tutorial, though these are not specific to meta-analysis.
2. with random effects in scale component
In the model with the random effects in the scale component, we can include random effects in the scale component, allowing residual heterogeneity to vary across studies. If we want to specify such a model, we can do so as follows:
sigma ~ 1 + Org.fertilizer.type + (1 | study_ID)This specification allows the residual variance to vary across studies, in addition to allowing study-level variation in the mean effect size via the location component. Conceptually, the location component models how the mean effect size varies with the moderator and across studies. The scale component models how the amount of unexplained variability differs across studies and fertiliser types.
Random effects in the location and scale components are assumed to be independent in this formulation.
3. Correlated random effects in location and scale components
The third model further extends the previous formulation by allowing correlated random effects between the location and scale components. For example, if we want to consider potential correlations between study-level deviations in mean effect sizes and deviations in residual variance, we can specify the model as follows:
(1 | p | study_ID)This specification allows the model to estimate the correlation between study-level deviations in the mean effect size and deviations in residual variance.
For example, a negative correlation would indicate that studies with larger mean effects tend to show lower residual variability, whereas a positive correlation would suggest the opposite. Such relationships are conceptually appealing and biologically interpretable in some contexts.
However, this model is computationally demanding and highly sensitive to prior specification and data structure. In applied meta-analyses, it frequently exhibits convergence problems or weakly identified correlation parameters.
glmmTMB
In contrast to metafor, glmmTMB allows us to include random effects in both location and scale components - we demonstrate how to fit a location-scale meta-regression model using glmmTMB.
But it should be mentioned that glmmTMB does not support correlated random effects between location and scale components, so the third model presented above cannot be implemented using glmmTMB.
Again, the same as brms, we focus on a location–scale meta-regression model with fixed effects in the scale component.
A. Only fixed effects in scale component
In glmmTMB, the location–scale structure is implemented by separating the model into a conditional component and a dispersion component. The conditional model specifies the mean structure of the response (location), while the dispersion model specifies how the residual variance changes as a function of predictors (scale). Note that the dispersion model in glmmTMB does not describe an additional random effect, but rather models heteroscedasticity in the residual variance on the log scale.
As shown below, the location–scale models with glmmTMB produce results that match those from brms completely.
lnRR
glmm_lnRR1 <- glmmTMB(lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
dispformula = ~ 1 + Org.fertilizer.type, # scale component
data = dat,
REML = TRUE)
summary(glmm_lnRR1)
# Family: gaussian ( identity )
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Dispersion: ~1 + Org.fertilizer.type
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 230.9 249.7 -110.5 220.9 313
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.0608129 0.24660
# g effect_size_ID33 0.0047213 0.06871
# effect_size_ID34 0.0426425 0.20650 0.00
# effect_size_ID35 0.0698313 0.26426 0.00 0.00
# effect_size_ID36 0.0798567 0.28259 0.00 0.00 0.00
# effect_size_ID37 0.0837435 0.28938 0.00 0.00 0.00 0.00
# effect_size_ID39 0.0095870 0.09791 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.0091520 0.09567 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.0363323 0.19061 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.0331452 0.18206 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.0012338 0.03513 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual NA NA
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.26791 0.05379 -4.981 6.33e-07 ***
# Org.fertilizer.typeplant -0.13991 0.05897 -2.373 0.0177 *
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Dispersion model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -1.35445 0.07922 -17.097 <2e-16 ***
# Org.fertilizer.typeplant 0.17055 0.09928 1.718 0.0858 . lnCVR
library(glmmTMB)
glmm_lnCVR1 <- glmmTMB(lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
dispformula = ~ 1 + Org.fertilizer.type,
data = dat,
REML = TRUE)
summary(glmm_lnCVR1)
# Family: gaussian ( identity )
# Formula: lnCVR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 + effect_size_ID | g, vcv)
# Dispersion: ~1 + Org.fertilizer.type
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 930.3 949.1 -460.2 920.3 313
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.1083168 0.32912
# g effect_size_ID33 0.0047213 0.06871
# effect_size_ID34 0.0426425 0.20650 0.00
# effect_size_ID35 0.0698313 0.26426 0.00 0.00
# effect_size_ID36 0.0798567 0.28259 0.00 0.00 0.00
# effect_size_ID37 0.0837435 0.28938 0.00 0.00 0.00 0.00
# effect_size_ID39 0.0095870 0.09791 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.0091520 0.09567 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.0363323 0.19061 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.0331452 0.18206 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.0012338 0.03513 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual NA NA
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 0.32787 0.09283 3.532 0.000413 ***
# Org.fertilizer.typeplant 0.07193 0.15208 0.473 0.636203
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Dispersion model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.41563 0.07030 -5.913 3.37e-09 ***
# Org.fertilizer.typeplant 0.66519 0.08807 7.553 4.26e-14 ***B. Random effects in location and scale components
If you want to include random effects in the scale component, you just need to insert them to the dispformula part of the model.
glmm_lnRR2 <- glmmTMB(lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) +
equalto(0 + effect_size_ID|g, vcv),
dispformula = ~ 1 + Org.fertilizer.type + (1 | study_ID),
data = dat,
REML = TRUE)
summary(glmm_lnRR2)
# Family: gaussian ( identity )
# Formula: lnRR ~ 1 + Org.fertilizer.type + (1 | study_ID) + equalto(0 +
# effect_size_ID | g, vcv)
# Dispersion: ~1 + Org.fertilizer.type + (1 | study_ID)
# Data: dat
#
# AIC BIC logLik -2*log(L) df.resid
# 177.4 199.9 -82.7 165.4 312
#
# Random effects:
#
# Conditional model:
# Groups Name Variance Std.Dev. Corr
# study_ID (Intercept) 0.0570683 0.23889
# g effect_size_ID33 0.0047213 0.06871
# effect_size_ID34 0.0426425 0.20650 0.00
# effect_size_ID35 0.0698313 0.26426 0.00 0.00
# effect_size_ID36 0.0798567 0.28259 0.00 0.00 0.00
# effect_size_ID37 0.0837435 0.28938 0.00 0.00 0.00 0.00
# effect_size_ID39 0.0095870 0.09791 0.00 0.00 0.00 0.00 0.00
# effect_size_ID40 0.0091520 0.09567 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID42 0.0363323 0.19061 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID44 0.0331452 0.18206 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# effect_size_ID45 0.0012338 0.03513 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
# ... ... ... ... ... ... ... ... ... ... ... ... ...
# Residual NA NA
# Number of obs: 318, groups: study_ID, 36; g, 1
#
# Dispersion model:
# Groups Name Variance Std.Dev.
# study_ID (Intercept) 0.06081 0.2466
# Residual NA NA
# Number of obs: 318, groups: study_ID, 36
#
# Conditional model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.26476 0.04887 -5.417 6.05e-08 ***
# Org.fertilizer.typeplant -0.11847 0.04545 -2.607 0.00914 **
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Dispersion model:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -1.6225 0.1147 -14.14 < 2e-16 ***
# Org.fertilizer.typeplant 0.4022 0.1406 2.86 0.00424 ** Software and package versions
R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.7.3
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Edmonton
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] patchwork_1.3.2 here_1.0.2 orchaRd_2.1.3
[4] tidybayes_3.0.7 lubridate_1.9.4 forcats_1.0.1
[7] stringr_1.6.0 dplyr_1.1.4 purrr_1.2.0
[10] readr_2.1.6 tidyr_1.3.2 tibble_3.3.0
[13] ggplot2_4.0.1 tidyverse_2.0.0 metafor_4.8-0
[16] numDeriv_2016.8-1.1 metadat_1.4-0 Matrix_1.7-4
[19] glmmTMB_1.1.14 brms_2.23.0 Rcpp_1.1.1
loaded via a namespace (and not attached):
[1] Rdpack_2.6.4 gridExtra_2.3 inline_0.3.21
[4] phangorn_2.12.1 sandwich_3.1-1 rlang_1.1.7
[7] magrittr_2.0.4 otel_0.2.0 matrixStats_1.5.0
[10] compiler_4.5.2 mgcv_1.9-3 loo_2.9.0
[13] vctrs_0.6.5 quadprog_1.5-8 pkgconfig_2.0.3
[16] arrayhelpers_1.1-0 fastmap_1.2.0 backports_1.5.0
[19] labeling_0.4.3 rmarkdown_2.30 tzdb_0.5.0
[22] nloptr_2.2.1 ggbeeswarm_0.7.3 xfun_0.55
[25] jsonlite_2.0.0 parallel_4.5.2 R6_2.6.1
[28] stringi_1.8.7 RColorBrewer_1.1-3 StanHeaders_2.32.10
[31] boot_1.3-32 estimability_1.5.1 rstan_2.32.7
[34] knitr_1.51 zoo_1.8-15 pacman_0.5.1
[37] bayesplot_1.15.0 splines_4.5.2 igraph_2.2.1
[40] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.17.1
[43] abind_1.4-8 yaml_2.3.12 TMB_1.9.19
[46] codetools_0.2-20 pkgbuild_1.4.8 lattice_0.22-7
[49] withr_3.0.2 bridgesampling_1.2-1 S7_0.2.1
[52] posterior_1.6.1 coda_0.19-4.1 evaluate_1.0.5
[55] RcppParallel_5.1.11-1 ggdist_3.3.3 pillar_1.11.1
[58] tensorA_0.36.2.1 checkmate_2.3.3 stats4_4.5.2
[61] reformulas_0.4.3.1 distributional_0.5.0 generics_0.1.4
[64] rprojroot_2.1.1 mathjaxr_2.0-0 hms_1.1.4
[67] rstantools_2.5.0 scales_1.4.0 minqa_1.2.8
[70] xtable_1.8-4 glue_1.8.0 emmeans_2.0.1
[73] tools_4.5.2 lme4_1.1-38 mvtnorm_1.3-3
[76] fastmatch_1.1-6 grid_4.5.2 ape_5.8-1
[79] rbibutils_2.4 QuickJSR_1.8.1 colorspace_2.1-2
[82] nlme_3.1-168 beeswarm_0.4.0 vipor_0.4.7
[85] cli_3.6.5 latex2exp_0.9.6 svUnit_1.0.8
[88] Brobdingnag_1.2-9 gtable_0.3.6 digest_0.6.39
[91] htmlwidgets_1.6.4 farver_2.1.2 htmltools_0.5.9
[94] lifecycle_1.0.5 MASS_7.3-65 References
- Lauren C. Ponisio, Leithen K. M’Gonigle, Kevi C. Mace, Jenny Palomino, Perry de Valpine, Claire Kremen; Diversification practices reduce organic to conventional yield gap. Proc Biol Sci. 2015; 282 (1799): 20141396. https://doi.org/10.1098/rspb.2014.1396